Under the risk neutral probability measure $\mathbb{Q}$, the underlying is defined as a geometric Brownian motion solving the following SDE:
\begin{equation} \begin{cases} \frac{d S_t}{S_t} = r dt + \sigma d W_t \\ S_0 = x \end{cases} \end{equation}From that SDE, one can explcitiely compute the first two moments of $S$:
\begin{equation} \begin{cases} \mathbb{E}_t^{\mathbb{Q}} (S_{t+h}) = S_t e^{r h} \\ \mathbb{E}_t^{\mathbb{Q}} (S_{t+h}^2) = S_t^2 e^{2 r h + \sigma^2 h} \end{cases} \end{equation}In the sequel, we will use a recombining tree to model the underlying moves. We define by $u, u > 1$, some multiplicative factor of the underlying from time $t$ to time $t+h$. In the trinomial model, the underlying at time $t+h$ can take 3 values knowing the $S_t$:
\begin{equation} \begin{cases} S_{t+h} = u \times S_t > S_t \text { with probability } p_u \\ S_{t+h} = d \times S_t < S_t \text { with probability } p_d \\ S_{t+h} = S_t \text { with probability } p_m = 1 - p_u - p_d \\ \end{cases} \end{equation}with by vertue of our recombining assumption, $d = \frac{1}{u}$.
That is, given $u$, we must find the matching probabilities $p_u, p_d, p_m$ matching the Black-Scholes model. Thus we have 3 unknown variables. We must at least give 3 independent equations to solve the system.
The probabilities should be non negative and summing to one, which is already written for the definition of $p_m$. We are left with 2 unknown variables $p_u$ and $p_d$. Matching the first two moments given above we have:
\begin{equation} \begin{cases} \mathbb{E}_t^{\mathbb{Q}} (S_{t+h}) = S_t e^{r h} = S_t \left(u p_u + (1 - p_u - p_m) + d p_d \right) \\ \mathbb{E}_t^{\mathbb{Q}} (S_{t+h}^2) = S_t^2 e^{2 r h + \sigma^2 h} = S_t^2 \left(u^2 p_u + (1 - p_u - p_m) + d^2 p_d \right) \end{cases} \end{equation}So we have:
$$ \begin{pmatrix} e^{r h} - 1 \\ e^{2rh + \sigma^2 h} - 1 \end{pmatrix} = \begin{pmatrix} u - 1 & d - 1 \\ u^2 - 1 & d^2 - 1 \end{pmatrix} \begin{pmatrix} p_u \\ p_d \end{pmatrix} $$or equivalently, solving the linear system:
And using the classical formula of the inverse of a 2x2 matrix, we have: $$ \begin{pmatrix} u - 1 & d - 1 \\ u^2 - 1 & d^2 - 1 \end{pmatrix}^{-1} = \frac{1}{(u - 1) (d^2 - 1) - (u^2 - 1) (d - 1)} \begin{pmatrix} d^2 - 1 & 1 - d \\ 1 - u^2 & u - 1 \end{pmatrix} $$
and $p_m := 1 - p_u - p_d$. Moreover, once the computation is done, one should check that all the probabilities are non negative.
In [1]:
import numpy as np
from scipy.stats import norm
In [2]:
class TrinomialBSModel(object):
def __init__(self, S0=100., r=0.02, sigma=0.2, mat=1.):
self.__s0 = S0
self.__r = r
self.__sigma = sigma
self.__T = mat
def __compute_probs(self):
B = np.array([-1. + np.exp(self.__r * self.__h),
-1. + np.exp(2. * self.__r * self.__h + self.__sigma**2 * self.__h)])
d = self.__down
u = self.__up
A = np.array([[u - 1., d - 1.],
[u**2 - 1., d**2 - 1.]])
det = (u - 1.) * (d**2 - 1.) - (u**2 - 1.) * (d - 1.)
invA = 1. / det * np.array([[d**2 - 1., 1. - d],
[1. - u**2, u - 1.]])
res = invA.dot(B)
self.__pu = res[0]
self.__pd = res[1]
self.__pm = 1. - self.__pu - self.__pd
assert 0 <= self.__pu <= 1., 'p_u should lie in [0, 1] given %s' % self.__pu
assert 0 <= self.__pd <= 1., 'p_d should lie in [0, 1] given %s' % self.__pd
assert 0 <= self.__pm <= 1., 'p_m should lie in [0, 1] given %s' % self.__pm
def __check_up_value(self, up):
if up is None:
lbda = np.sqrt(0.5 * np.pi)
up = np.exp(lbda * self.__sigma * np.sqrt(self.__h))
assert up > 0., 'up should be non negative'
down = 1. / up
assert down < up, 'up <= 1. / up = down'
self.__up = up
self.__down = down
def __gen_stock_vec(self, nb):
vec_u = self.__up * np.ones(nb)
np.cumprod(vec_u, out=vec_u)
vec_d = self.__down * np.ones(nb)
np.cumprod(vec_d, out=vec_d)
res = np.concatenate((vec_d[::-1], [1.], vec_u))
res *= self.__s0
return res
def payoff(self, stock_vec):
raise NotImplementedError()
def compute_current_price(self, crt_vec_stock, nxt_vec_prices):
expectation = np.zeros(crt_vec_stock.size)
for i in range(expectation.size):
tmp = nxt_vec_prices[i] * self.__pd
tmp += nxt_vec_prices[i + 1] * self.__pm
tmp += nxt_vec_prices[i + 2] * self.__pu
expectation[i] = tmp
return self.__discount * expectation
def price(self, nb_steps, up=None):
assert nb_steps > 0, 'nb_steps shoud be > 0'
nb_steps = int(nb_steps)
self.__h = self.__T / nb_steps
self.__check_up_value(up)
self.__compute_probs()
self.__discount = np.exp(-self.__r * self.__h)
final_vec_stock = self.__gen_stock_vec(nb_steps)
final_payoff = self.payoff(final_vec_stock)
nxt_vec_prices = final_payoff
for i in range(1, nb_steps + 1):
vec_stock = self.__gen_stock_vec(nb_steps - i)
nxt_vec_prices = self.compute_current_price(vec_stock, nxt_vec_prices)
return nxt_vec_prices[0]
In [3]:
class TrinomialBSCall(TrinomialBSModel):
def __init__(self, S0=100., r=0.02, sigma=0.2, mat=1., K=100.):
super(TrinomialBSCall, self).__init__(S0, r, sigma, mat)
self.__K = K
def payoff(self, s):
return np.maximum(s - self.__K, 0.)
class TrinomialBSAmericanCall(TrinomialBSCall):
def compute_current_price(self, crt_vec_stock, nxt_vec_prices):
crt_payoff = self.payoff(crt_vec_stock)
crt_prices = super(TrinomialBSAmericanCall, self).compute_current_price(crt_vec_stock, nxt_vec_prices)
return np.maximum(crt_payoff, crt_prices)
In [4]:
class TrinomialBSPut(TrinomialBSModel):
def __init__(self, S0=100., r=0.02, sigma=0.2, mat=1., K=100.):
super(TrinomialBSPut, self).__init__(S0, r, sigma, mat)
self.__K = K
def payoff(self, s):
return np.maximum(self.__K - s, 0.)
class TrinomialBSAmericanPut(TrinomialBSPut):
def compute_current_price(self, crt_vec_stock, nxt_vec_prices):
crt_payoff = self.payoff(crt_vec_stock)
crt_prices = super(TrinomialBSAmericanPut, self).compute_current_price(crt_vec_stock, nxt_vec_prices)
return np.maximum(crt_payoff, crt_prices)
In [5]:
def bs_call_price(S=100., r=0.02, sigma=0.2, t=0., T=1., K=100.):
ttm = T - t
if ttm < 0:
return 0.
elif ttm == 0.:
return np.maximum(S - K, 0.)
vol = sigma * np.sqrt(ttm)
d_minus = np.log(S / K) + (r - 0.5 * sigma**2) * ttm
d_minus /= vol
d_plus = d_minus + vol
res = S * norm.cdf(d_plus)
res -= K * np.exp(-r * ttm) * norm.cdf(d_minus)
return res
In [6]:
def bs_put_price(S=100., r=0.02, sigma=0.2, t=0., T=1., K=100.):
# Using call-put parity :)
ttm = T - t
if ttm < 0:
return 0.
elif ttm == 0.:
return np.maximum(K - S, 0.)
dsct_strike = K * np.exp(-r * ttm)
cap_S = S * np.exp(r * t)
call = bs_call_price(S, r, sigma, t, T, K)
return call - cap_S + dsct_strike
In [7]:
tree = TrinomialBSCall()
print(tree.price(1000))
print(bs_call_price())
In [8]:
tree = TrinomialBSPut()
print(tree.price(1000))
print(bs_put_price())
In [9]:
tree = TrinomialBSAmericanCall()
print(tree.price(1000))
In [10]:
tree = TrinomialBSAmericanPut()
print(tree.price(1000))